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Abstract. In this paper we present an algorithm for pricing barrier options in one-dimensional 
Markov models. The approach rests on the construction of an approximating continuous-time 
Markov chain that closely follows the dynamics of the given Markov model. We illustrate the 
, method by implementing it for a range of models, including a local Levy process and a local 

' volatility jump-diffusion. We also provide a convergence proof for this algorithm. 

(N 

^ : 

' 1. Introduction 

l> . 

I 1.1. Background and motivation. Barrier option are among the most popular exotic derivatives. 
Such contracts form effective risk management tools, and are very liquid in the Foreign Exchange 

Pm ■ 

p ^ , markets (see e.g. Wystup [5D], Hakala and Wystup P^). The most liquid barrier options in FX 
, markets are continuously monitored single- or double-no-touch options and knock-in or knock-out 

■ calls and puts (see e.g. Lipton [40j . |41j). The main challange in the risk management of large 
^ portfolios of barrier options faced by trading desks that make a market in these securities is to 

be able to price and hedge the barrier products in models that are flexible enough to describe the 

^ ' observed option prices (i.e. calibrate to the vanilla market). 
00 ' 

' It is by now well established that the classical Black-Scholes model lacks the flexibility accurately 

I to fit to observed option price data (see e.g. Gatheral [20] and the references therein). A variety 

0^ I of models have been proposed to provide an improved description of the dynamics of the price of 
i the underlying that can more accurately describe the option surface. Parametric diffusion models 
, like the CEV process [13] have additional flexibility to fit the vanilla skew at a single maturity for 
J> , as many options as there are free parameters in the model. The seminal idea (see Dupire |17j and 

^ ■ Gyongy [22]) that allows one to construct a model that can describe the entire implied volatility 

■ surface (across all strikes and maturities) is that of local volatility models, where a non-parametric 
form of the local volatility function is constructed from the option price data. It has been shown 
that in practice such models imply unrealistic dynamics of the option prices (see the formula for 
the implied volatility in a local volatility model given in [23]). The ramification is an unrealistic 
amount of vega risk, which is very expensive to hedge. Therefore, even though in a local volatility 
model barrier options can be priced using a PDE solver, this modelling framework alone is not 
suitable for the risk management of a large portfolio of barrier options. 

At the other end of the spectrum are the jump processes with stationary and independent 
increments, which can fit very well the volatility smile at a single maturity (see e.g. [12] and the 
references therein). A variety of models in the exponential Levy class have been proposed in the 
literature: CGMY [8], KoBoL [6], generahsed hyperbolic [IB], NIG [1] and Kou [32]. Exponential 
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Levy processes are simple examples of Markov processes whose law is uniquely determined by 
the distribution of the process at a fixed time. Since the set of call option prices at a fixed 
maturity for all strikes uniquely determines the marginal risk-neutral distribution at that maturity, 
calibration to option prices at multiple maturities in principle fixes the corresponding marginals. It 
has been observed (see e.g. [12j) that Levy processes lack the flexibility of calibrating simultaneously 
across a range of strikes and maturities. Several generalisations within the one-dimensional Markov 
framework have been proposed. 

If the stationarity assumption is relaxed while the independence of increments is retained, one 
arrives at the class of exponential additive processes, which have recently been shown to calibrate 
well to several maturities in equity markets. The Sato process introduced in Carr et al. [TO] is an 
example of such an additive model used in financial modelling. 

The independent increments property of a process implies that its transition probabilities are 
translation invariant in space, so that they only depend on the difference between the end and 
starting value of the process. It is well known that the distribution of a log-asset price depends 
in a non-linear way on the starting point (e.g. in equity markets it has been observed that if the 
current price is high, then the volatility is low and vice versa). To capture this effect one is led 
to consider Markov jump-processes whose increments are not independent. As a generalisation of 
local volatility models, the class of local Levy processes introduced by Carr et al. [S] allows the 
modeller to modulate the intensity of the jumps as well as their distribution depending on where the 
underlying asset is trading. A local volatility jump-diffusion with similar structural properties was 
calibrated to the implied volatility surface in Andersen and Andreasen [3] and He et al. [25j). Due to 
the presence of jumps and the absence of stationarity and independence of increments, the problem 
of obtaining the first-passage probabilities for such a general class of processes is computationally 
less tractable. 

There exists currently a good deal of literature on numerical methods for the pricing of barrier- 
type options. It is well known that in this case a straightforward Monte Carlo simulation algorithm 
will be time-consuming and yield unstable results for the prices and especially the sensitivities. The 
knock- in/out features in the barrier-option payoffs also lead to slower convergence of the Monte 
Carlo algorithm. To address this problem the following (semi-)analytical approaches have been 
developed for specific models: 

(a) spectral expansions for several parametric diffusion models (Davydov and Linetsky [14] . 
Lipton [40J), 

(b) transform based approaches for exponential Levy models (Boyarchenko and Levendorskii [3], 
Geman and Yor [21], Jeannin and Pistorius [30], Kou and Wang [33j, Sepp [49] ) . 

The method (a) employs the explicit spectral decompositions for this class of diffusion models, 
whereas the approach (b) exploits the independence and stationarity of the increments of the Levy 
process, and the so-called Wiener-Hopf factorisation. Since both of these approaches hinge on 
special structural properties of the underlying processes, it is not clear how (and if) they can be 
extended to more general Markovian models. 
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Figure 1. This is a schematic picture in bloc notation that demonstrates how the matrices A 
and A are obtained from the generator A. The matrix A is the generator matrix of an approximating 
continuous-time Markov chain X on the state-space G. The subset G C G consists of the elements 
of G that lie between the barriers. The matrix A contains all the information necessary to price 
any contract that knocks out (or knocks in) when a barrier is breached. Similarly the matrix A is 
what is required to compute the distribution of any function that depends on the first-passage and 
overshoot of the chain X into the region on the other side of the barriers. 



A different approach, pioneered by Kusliner (see e.g. [36]), is the discrete time Markov chain 
approximation method. Originally developed for the numerical solution of stochastic optimal con- 
trol problems in continuous time, this method consists of approximating the system of interest by 
a chain that closely follows its dynamics and solving the problem of interest for the chain. An 
application to the pricing of American type options was given in [35]. Duan et al. [15] priced a 
discretely monitored barrier option in the Black-Scholes and NGARCH models, using a discrete 
time Markov chain. Rogers and Stapleton [36] develop an efficient binomial tree method for barrier 
option pricing (see also references therein for related methods). Zvan et al. [51] investigate a PDE 
finite difference discretization methods for barrier and related options. Markov chains have also 
been employed as a modelling tool for price processes; Albanese and Mijatovic [2J modelled the 
stochasticity of risk reversals and carried out a calibration study in FX markets under a certain 
continuous-time Markov chain constructed to model the FX spot process. 

1.2. Contribution of the current paper. In this paper we consider the problem of pricing 
barrier option in the setting of one-dimensional Markov processes, which in particular includes 
the case of local Levy models as well as local volatility jump-diffusions and additive processes. 
The presented approach is probabilistic in nature and is based on the following two elementary 
observations: (i) given a Markov asset price process S it is straightforward to construct a continuous- 
time Markov chain X whose law is close to that of S, by approximating the generator of the process 
S with an intensity matrix; (ii) the corresponding first-passage problem for a continuous-time 
Markov chain can be solved explicitly via a closed-form formula that only involves the generator 
matrix of the chain X. More precisely, for a given Markov asset price model S on the state-space 
E = (0, cxd) with corresponding generator C the algorithm for the pricing of any barrier product 
(including rebate options, which depend on the position at the moment of first-passage) consists of 
the following two steps. 
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(1) Construct a finite state-space G C E and a generator matrix A for the chain X that 
approximates the operator C on G. 

(2) To value knock-out and rebate options, obtain the matrices A and A by the procedure in 
Figure [1] and apply the closed- form formulas in equations (j3.7p and (j3.8p . 

The form of the generators of Markov processes that commonly arise in pricing theory (including 
the local Levy class) is well known from general theory (they are reviewed in Sections [2] and H]) . 

The state-space G in step (1) is generated using a standard procedure from the PDE literature 
(described in Appendix [B]). The generator matrix A is defined by matching the instantaneous local 
moments of the Markov processes S and the chain X on the state-space G. This criterion implies 
in particular that the chain X locally drifts at the same rate as the asset price process S. 

Step (2) of the algorithm consists of the evaluation of the closed-form formulas for the first- 
passage probabilities that can be derived employing continuous-time Markov chain theory (Theo- 
rem[T]). The evaluation of this formula consists of exponentiation of either matrix A or A, which can 
be performed using the Fade approximation algorithm that is implemented in standard packages 
such as Matlab (see also [26]). 

We implemented this algorithm for a number of models, and found good agreement with the 
numerical results obtained for models considered elsewhere in the literature (see Section [7j) . We 
also prove that by refining the grid the prices generated by this approach converge to those of 
the limiting model (see Section ^ . There is much literature devoted to the study of the (weak) 
convergence of Markov chains to limiting processes. However the (exact) rates of convergence of 
prices generated by Markov chains to those of the limiting model are rarely available, especially for 
barrier options. Establishing the rates for specific models remains an open question, left for future 
research. 

The remainder of the paper is organized as follows. In Section [2] we define the precise class of 
models and barrier option contracts that is considered, and state some preliminary results about 
Markov processes. Section [3] presents the formulas for the first-passage quantities of the continuous- 
time Markov chains. In Section[3]we extend the approach to time-inhomogeneous Markov processes. 
In Section [5] we describe the discretization algorithm that constructs the intensity matrix A of the 
chain X. Section [6] states the convergence results, which are proved in Appendix lAl Numerical 
results are presented in Section [7] and Section [8] concludes the paper. 

Acknowledgements: we would like to thank Dilip Madan and the participants of the 2009 Leicester 
workshop on Spectral and Cubature Methods in Finance and Econometrics for useful suggestions. 
Research supported by EPSRC grant EP/D039053. 

2. Problem setting: Barrier options for Markov processes 

The problem under consideration is that of the valuation of general barrier options, which can 
be formulated as follows. Given a random process S = {St}t>o modelling the price evolution of 
a risky asset, non- negative payoff and rebate functions g and /i, and a set A specifying the range 
of values for which the contract 'knocks out', it is of interest to evaluate the expected discounted 
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value of the random cash flow associated to a general barrier option contract 

(2-1) g{ST)I{r^>T} + hiSrJI{r^<T}, 

where Ic denotes the indicator of a set C and 

TA = inf{t >0: St£ A} 

is the first time that S enters the set A. Furthermore, it is relevant to quantify the sensitivities 
of this value with respect to different parameters such as the spot value = x. The cash flow 
in ()2.ip consists of a payment g^Sx) in the case the contract has not knocked out by the time T, 
and a rebate h(Srj^ ) if it has. Examples of commonly traded options included in this setting are 
the down-and-out, up-and-out and double knock-out options. In particular, by taking j4 = we 
retrieve the case of a standard European claim with payoff h{ST) at maturity T. 

We will consider this valuation problem in a Markovian setting, assuming that the underlying 
is a Markov process with state-space E := (0,oo) defined on some filtered probability space 
F,P), where F = {J^t}t>o denotes the standard filtration generated by S. Thus, S takes 
values in E and satisfies the Markov property: 

(2.2) E[fiSt+s)\:Ft]=Psf{St), 

for all s,t > and bounded Borel functions /, where E denotes the expectation under the proba- 
bility measure P and Psf is given by 

(2.3) Psfix) := E,.[/(5,)] := E[/(5,)|5o = x]. 

By taking expectations in (j2.2p we see that the family {Pt)t>o forms a semgroup: 

Pt+sf = PtiPsf), for all s,t>0, and Pq/ = /• 

Informally, these conditions state that the expected value of the random cash flow f{St+s) occurring 
at time t + s conditional on the available information up to time t depends on the past via the 
value St only. Setting the rate of discounting equal to a non-negative constant r, for any pair of 
non-negative Borel functions g and h the expected discounted value of the barrier cash flow (j2.1|) 
at the epoch ta A T, the earlier of maturity T and the first entrance time ta, is given by 

(2.4) [e-^^g{ST)I{r^>T}] + [e— ^/i(5.JI{.^<t}] . 

If S represents the price of a tradeable asset, r is the risk-free rate, d is the dividend yield and the 
process {e~^^~'^^^ St}t>o is a martingale, standard arbitrage arguments imply that no arbitrage is 
introduced if expression ()2.4p is used as the current price of the payoff ()2.ip . 

Before proceeding we will review some key concepts of the standard Markovian setup that will 
be needed. For background on the (general) theory of Markov processes we refere to the classical 
works Chung and Walsh [TT], Ethier and Kurtz 19], Ito and McKean [28j and Rogers and Williams 
[47j (the latter two in particular treat Markov processes in the context of diffusion theory). In 
what follows we will restrict S to be in a subclass of Markov processes for which, if the function / 
is continuous and tends to zero at infinity, the expected payoff Ptf{x) has the following properties: 
it depends continuously on the present value Sq = x and on expiry t and also decays to zero when 
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X tends to infinity. More precisely, denoting by Co(E) the set of continuous functions / on E that 
tend to zero at infinity and at zero, we make the following assumption: 

Assumption 1. S is a Feller process on E, that is, for any f € Co(E), the family {Ptf)t>Q, with 
Pff defined in (12. 3p . satisfies: 

(i) Ptf e Co(E) for any t > 0; 

(ii) limj^o Ptf{x) = f{x) for any x G E. 

The Feller property is a standard condition, which guarantees that a version of the process S 
with cadlag paths exist, satisfying the strong Markov property. In particular, a Feller process is a 
Hunt process. 

Throughout the paper we will take the knock-out set A to be of the form 
(2.5) A = (0, £] U [u, oo), < ^ < M < oo, 

which includes the cases of double and single barrier options — the latter by taking £ = or li = oo. 
To rule out degeneracies we will make the following assumption on the behaviour of S at the 
boundary points I and u: 

Assumption 2. 'Px{ta = t^o) = 1 where A° = {0,£) U (n, oo). 

This assumption states that the first entrance times into A and its interior coincide. A sufficient 
condition for Assumption [2] to be satisfied is Pa;(r^o = 0) = 1 for x G {£,u}', that is, when 
started at i or n, the process S immediately enters the interior of A. The class of Feller processes 
satisfying Assumption [2] includes many of the models employed in quantitative finance such as 
(Feller-) diffusions, jump-diffusions and Levy processes whose Levy measure admits a density. 

The family {Pt)t>o is determined by its infinitesimal generator C which is a map defined in the 
following way. Let P C Co(E) be the set of all / G Co(E) for which the right-hand side of (j2.6|) 
converges in the strong senseo Then T) is dense in Cq(E) and for any / G P, the function Cf is 
defined as 



These fundamental facts about semigroups and their generators can be found in [19\ Ch. 1]. 
We next give a few examples of Feller processes with their generators. 

Example 1. In a diffusion the asset price model S = {St}t>o evolves under a risk- neutral measure 
according to the stochastic differential equation (SDE) 



where So G E is the initial price, 7 G M and a : M+ R-^- is a given measurable function. To 
guarantee the absence of arbitrage we assume that a is chosen such that the discounted process 
{e~'^^St}t>o is a martingale, which implies in particular that S does not explode to infinity. If, in 

"'^That is, the convergence is with respect to the norm ||/|| := sup^g^ \f{^)\ oi the Banach space (Co(E), || • ||). 



(2.6) 



£f{x) :=liml(Pi/-/)(x). 



(2.7) 



—l = jdt + a {St) dWt 
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addition, infinity is not entrancqj for S and a is a continuous function, then 5 is a Feller process, 
and its infinitesimal generator acts on / € C^(Ej§ as 

^/(^) = ^^^(A/)(x) +7x(V/)(x), 
where V/ and A/ denote the first and second derivatives of / with respect to x (see [121 Sec.8.1]). 
Example 2. (a) The price process in an exponential Levy model S given by 



e 



Lt 



where r and d are constants representing the interest rate and dividend yield and L = {Lt}t>o 
is a Levy process, such that Eo[e^*] < oo for all t > 0. By construction, {e~*^^~'^)*S't}f>o is a 
martingale. The law of L is determined by its characteristic exponent ^, which is defined in terms 
of the characteristic function of Lt by = exp(t^'(u)) and, according to the Levy-Khintchine 

representation, takes the form 

^{u) = leu - -— + / (e^"2^ - 1 - u{dy), 

where {c,a'^,i') is the characteristic triplet, with o", c G M and u the Levy measure, which satisfies 
the integrability condition J^{1 A \y\'^)i'{dy) < oo. Further, Eo[e^*] < oo if and only if the Levy 
measure v integrates exp(y) at infinity, that is, 

/oo 
e^uidy) < oo. 

The process is a Feller process with an infinitesimal generator acting on / G C'c (^) i^^- Sato 
pi Thm. 3L5]) 



where 



^/(^) = ^A/(x) + CxVfix) + / ifixey) - fix) - xV/(x)(e^ - l)^y^^,}Hdy), 

^ J -oo 

/oo 
- l)I{|,|>i}Z.(dy). 
-oo 



(b) A closely related model is a geometric Levy process specified by the SDE 

—1 = (r - d - fi)dt + dLt, t > 0, 

where So,fi £ (0, oo) and L is a Levy process with characteristic triplet (c, ci^, z^), where the Levy 
measure ly has support in (— l,oo) and zvidz) < oo. The former condition guarantees that 
St > for all t > 0, the latter that E[|Lj|] < oo for all t > 0. The discounted process {e~^''~'^^^ St}t>o 
is a martingale if /i is given by 

/oo 
ziy{dz). 



See Ito and McKean for an explicit criterion in terms of 7 and a for this to be the case. 
'Cc(E) denotes the set of functions with compact support in E = (0, 00). 
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Further, 5 is a Feller process, with an infinitesimal generator L that acts on / G C^(E) as 
Cf{x) = -^Af{x) + CxVfix) + J ^ [f{x{l + z)) - fix) - xV/(x)zI|,<i|]i/(dz), 
where (^ = c + r — d — fi. For further background on Levy processes we refer to Sato 



Example 3. More generally, one may specify the price process S by directly prescribing its gen- 
erator C to act on sufficiently regular functions / as 

(2.10) Cfix) = ^^^Af{x) + {r-d-fi{x))xVf{x) 



/OO 
^ [/(x(l + y)) - fix) - V/(x)xyI|iy|<i}]z.(x, dy), 



where //,(T : E — > M are given functions, and for every x € E, i/(x,dy) is a (Levy) measure 
with support in (— l,oo) such that J^y^z^(x,dy) < oo. Sufficient conditions on fi,a and v to 
guarantee the existence of a Feller process S corresponding to this generator were established in 
Kolokoltsovy [3ll Thm. 1.1]. If it holds that 

/oo 
yvix,dy) < oo, 

then the discounted process St}t>o is a local martingale. 

The barrier option payoff (|2.1|) can be expressed in terms of the process S that is stopped when 
it enters the set A, which we denote by = StAr^ > as 

[e-'^^5(ST)I{r,>T}] + [e-'^-^/i(5.JI|.^<T}] = E,[e-'-(^^^)/(5TArJ] 

=: P^fix), 



where the function / is defined as 



/i(x), X £ A, 
gix), else. 



The process is itself a Markov process, and, as a consequence, the family iP^)t>o satisfies the 
semigroup property. The infinitesimal generator corresponding to can be explicitly expressed 
in terms of C as follows: 

Lemma 1. For any f £ D, where D is the domain of the generator C, we have 



(2.12) limt-i(Pt^/(x) - fix)) = kix) := { 



X e A 

(£-r)/(x) x(^A 



where the convergence is pointwise. If is itself a Feller process and the function f in addition 
satisfies lim^j-^a^ £/(x) = 0, where OA is the boundary of A, then C^f = k, where is the 
infinitesimal generator of the semigroup (P/^)j>o. 



CONTINUOUSLY MONITORED BARRIER OPTIONS UNDER MARKOV PROCESSES 9 

If is a Feller process, the relation between and £^ can formally be expressed as 
(2.13) = exp(f£^). 

Equation ()2.13p can be given a precise meaning if, for example, P/^ can be defined as a self-adjoint 
operator on a separable Hilbert space (see e.g. Ch. XII in Dunford and Schwarz [16], or Hille 
and Philips [23 )• determining the spectral decomposition of one can construct a spectral 
expansion of P/^f{x), which in the case of a discrete spectrum reduces to a series expansion. See 
Linetsky [67\ f551 f5U] for a development of this spectral expansions approach for one-dimensional 
diffusion models in finance, and an overview of related literature. 

When (asymmetric) jumps are present, the operator is non-local and not self-adjoint, and the 
spectral theory has been less well developed, and fewer explicit results are available. Here we will 
follow a different approach: we will approximate S by a finite state Markov chain, and show that 
for the approximating chain a matrix analog of the identities ()2.12p - (j2.13p holds true, where the 
infinitesimal generators C"^ can be easily obtained from £. In Section Owe give a self-contained 
development of this approach, and present an extension to time-dependent dynamics in Section [H 

3. Exit probabilities for continuous-time Markov chains 

Given a Markov price process S of interest, the idea is to construct a continuous-time finite 
state Markov chain X that is "close" to S, and to calculate the relevant expectations for this 
approximating chain. In this Section we will focus on the latter; we will return to the question of 
how to construct the chain in Section 6. Assume therefore we are given a continuous-time Markov 
chain X = {Xt}t>o- From Markov chain theory it is well known that the chain is completely 
specified by its state-space (or grid) G C E and its generator matrix A, which is an x square 
matrix with zero row sums and non-positive diagonal elements, if G has N elements. Given the 
generator matrix A the family of transition matrices {Pt)t>o of X, defined by Pt{x, y) := [Xt = y] 
for x,y £ G, is given by 

Pt=exp{tA). 

Conversely, the generator A can be retrieved from {Pt)t>o by differentiation at t = 0, that is, 

Pt = I + tA + o{t), 
as t I 0. Thus, for the chain X, the expected value of (j){Xt) is given by 
(3.1) E.MXt)] = (exp(tA)(/.) (x) 

for any function </> : G ^ M. Here and throughout the paper we will identify any square matrix 
A E M^^^ and any vector in with functions 

^ : G X G ^ M, A{x, y) := e'^Acy, x,y £ G, and 

: G ^ M, 4>{x) := e'^(t), x G G, 

where the vectors ex,ey denote the corresponding standard basis vectors of and ' stands for 
transposition. 
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Below we will show how to express the exit probabilities of the chain using matrix exponentiation 
in a way that is identical in form to the formula in ()3.ip . To that end, we partition G into a 
'continuation' set G and a 'knock-out' set := G\G, where 

(3.2) G 



and define the first exit time of X from G by 



(3.3) 



mi{t G 



where we use the convention inf := cx) and where we will take the set A as in ()2.5p . 

The value of a general barrier knock-out option with a rebate depends on the joint distribution of 
the exit time ta from A and the positions of the underlying at maturity and at the moment of exit. 
The corresponding quantities for the chain X can be expressed in terms of two transformations of 
X, X and X, namely the chain that is killed upon exiting G and the one that is absorbed at that 
instance, respectively. Correspondingly, we associate to the generator matrix A, two matrices: the 
N X N matrix A, where N := |G|, and the N x N matrix A^., defined by 



(3.4) 



(3.5) 



K{x,y) 



A{x,y) 



A{x,y) 
A(x,y) 


A{x,y) 



- r if X £ 

X £ 
X £ 
if X £ G, y £ 



, x = y, 
, y e G, X / y, 
^ y G G, 



where If; is the indicator function of the set G. We can now state the key result of this section: 



Theorem 1. For any T > 0, x £ G and r > and any function (p ■ 



it holds that 



(3.6) 



-r(TAT) 



0(X 



TAt 



exp ( TAr ] (j)] (x). 



In particular, for ip : G and ^ : G ^ M with ^(x) = for x £ G we have that 

(3.7) [4;{XT)I{r>T}] = (exp (ta) ?/^) (x) for any x £ G, 

(3.8) B^[e-^^aXr)Iir<T}] = (exp (ta.) (x) for any x £ G. 

Formulas (I3.7|) -( f3r8]l give a simple and robust way of computing barrier option prices by a 
single matrix exponentiation. The expectation in (13. 7p can be obtained by computing the spectral 
decomposition of the matrix A = UDU^^ and applying the formula exp (^A^ — U exp{TD)U^^ . 
The powerful Fade approximation method for matrix exponentiation, described in [26], can also be 
used to compute efficiently the matrix exponentials in Theorem [H particularly in the case where 
the matrices involved are of large dimension. Note that Theorem [T] can be seen to follow by an 
application of Lemma [H in order to clarify the ideas underlying the algorithm we present next a 
direct probabilistic derivation. 

Proof. To prove equation ()3.6p . we will verify that the expected value of an Arrow-Debreu barrier 
security that pays 1 precisely if X is in the state y at the earlier of the maturity T and the knock-out 
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time T is given by 

(3.9) E,[e-''(^^^)l{x^^,=j;}] = (exp(rA,)) {x, y) for all x,yeG. 

For a given time grid T„ = {kAt, k = 0,1,2, .. .) with At = T/n denote by Pt{x, y) the expected 
value of the corresponding discretely monitored Arrow-Debreu security and let 

Tn = inf{s eTn-.Xs^G} 

be the corresponding time at which the barrier is crossed. Since the paths of the chain X are 
piecewise constant, it follows that t„ | r and Xr,^ Xr as n tends to infinity. Hence the expected 
values of the discretely monitored Arrow-Debreu securities converge to the expected value of the 
continuously monitored one, 

PHx,y) = E..[e-(^--)l|^_^=,}] E.[e-'-(^-^)l|;,_=,j]. 

Clearly, since is the knock-out set, it holds for all i > that 

1 if X € G'^, X = y 
if X G G'^, X ^ y 
I-l){x,y) for all x G G^ y G G, 



Pt{x,y) 



where I is a square matrix of size with I[x, x) = 1 if x G G and zero else. Further, for x G G an 
application of the Markov property of X shows that 

PT{x,y) = e-^^*^PAt(2;,^)PT-At(2,y) 

lie-'-^'PAt) Pr-At) {x,y). 

Combining the two cases, iterating the argument and using the differentiability of Pt at t = shows 
that 

At , 



PT{x,y) = [[I-I + Ie-'^'P^t)PT-At){x,y) 
I-l + le-^^'PAtf^') {x,y) 



'l + l^-rAt^p^^ _ J) +7(e-At _ ^^-^T/At^ ^^^^^ 

/ ~ - \T/At\ 

{l + At(Ko-rI)+o{At)) Ux,y), 

since Aq = /A. When At = T/n tends to zero, this expression converges to ^exp ^TA^^^ {x,y), 
which completes the proof of (|3.9p and hence implies (j3.6p . Equation (j3.8p follows then directly by 
applying (j3.6p to ^. Finally, noting that for any -0 : G ^ M and any x G G we have 

{ktp^ (x) = (^AoV-o) {x), where i/'o : G ^ M is given by V'o(y) := '0(y)I{yge}, 
we get that 

(exp(rA)?/.) (x) = (exp(rAo)V'o) {x), 
which yields (13. 7p . |--| 
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4. TiME-INHOMOGENEOUS MARKOV PROCESSES 

In practical applications model parameters such as the short rate or the volatility function are 
often taken to be time-dependent. In this section we will briefly show how the results of the earlier 
section can be adapted to this time-dependent setting. 

We now start from a time-inhomogeneous Markov process S = {St}t<=[o,T]j modelling the evo- 
lution the risky asset price under consideration until some maturity T > 0. Then the time-space 
process {{t, St)}te[o,T] forms a two-dimensional homogeneous Markov process, and the approach 
developed in Section [2] can be easily adapted to this case, which we now outline. To highlight the 
role of time as a factor we denote hy D = {Dt}t>o and Y = {Yt}t>o where 

(4.1) Dt := {Do +t)AT and Yt := Sd^ Vt > 0. 

Then, for any s,t > and bounded Borel function / : [0, T] x E ^ R it holds that 

E[f{Dt,Yt)\J^s] = {Qtf)iDs,Ys), 

where 

{Qtf){s,x) = Bs,M{Dt,Yt)]:=B[f{Dt,Yt)\Do = s,Yo = x] 

= B[f{s + t,S^s+t)AT)\Ss = x]. 

Thus the family {Qt)t>o forms a (two-dimensional) semigroup. We will restrict ourselves to the 
case when (D, Y) is a Feller process with state-space E' := [0, T] x E, in the sense of Assumption [T] 
with E replaced by E'. Associated to the semigroup {Qt)t>o is the infinitesimal generator C' defined 

by 

\imt-'[Qtf{s,x) - f{s,x)] = £'f{s,x), 

for all / in the set T>' , the domain of for which this limit exists. The process stopped upon 
entering a closed set A remains a Markov process, and its corresponding generator can be obtained 
as in Section [2l Denote by Ct the generator restricted to functions 51 : E ^ M of space only, that 
is, {Ctg){x) := {C'g){t,x). 

We next give two examples of models in which time-dependence plays a natural role — we will 
present a numerical illustration of these models in Section [71 

Example 4. Given a surface of arbitrage-free call option prices generated by some diffusion model, 
there exists a one-dimensional diffusion with time-dependent volatility function E(t, x) that repro- 
duces those option prices. The volatility function S(t,x) is explicitly given in terms of the call 
option prices by Dupire's formula. In such a local volatility model the stock price S = {St}t&[Q,T] 
evolves according to the SDE 

(4.2) dSt = 7(t, St)dt + S(t, St)dWt, 

where 'y{t,x) = {r{t) — d{t))x with r{t) and d{t) continuous functions representing the interest 
and dividend rate, and Y,{t,x) is a volatility function chosen sufficiently regular that {{Dt,Yt)}t>o 
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is a Feller process and the discounted stock price {e~ St}te[o,T] is a martingale. The 

corresponding infinitesimal generator then acts on functions / G C'c(^) ^ 

{Ctf){x) = 7(t,x)V/(x) + is2(t,x)A/(x), 

where V/, A/ denote the first and second derivatives with respect to x. 

Example 5. In an exponential additive model the log-price is a spatially homogeneous Markov 
process, and S = {St}te[o,T] is given by 

gVt 

where X is an additive process, a process with independent increments and cadlag paths. If 
E[e''^*] < oo for all t G [0,T], the discounted process {e-i'o[^(")-'='(")1^^5t}te[o,T] is a martingale, by 
construction. 

To be specific, let X be an additive process X = {^t}tg[o,T] with drift and volatility functions 
P{t) and a{t), and a system of Levy measures dy), given by 



(4.3) St = e^oiHs)-d{s)]ds 



(4.4) 



Xt= [ (3{s)ds+ I a{s)dW, + Zt. 
Jo Jo 



Here the process Z = {Zt}te[o,T] is a jump-process with compensator i^{t,dy)dt = g{t,y)dydt 
satisfying the integrability condition 

(4.5) [ (lAy2)sup{5(s,y) :sG [0,r]}dy<cx), 

Jr\{o} 

with g{-,y) : [0, T] M+, y £ M\{0}, a family of continuous functions. Models from this class 
(with cr = and a particular form of g{t, y)) have been proposed for stock price modelling by Carr 
et al. [TO] — see also Section \77l\ Then the time-space process {{Dt,Yt)}t>o is a Feller process, 
with an infinitesimal generator acting on / G C|(E) aj^ 

(4.6) {Ctf){x) = (£7)(t,x)+ / [/(xe^)-/(x)-V/(x)x(e^-l)I{|j^l<i}]5(t,y)d2/, 



(£7)(t,x) := P{t)xVf{x) + ^x^Af{x), 

where 

m = r{t) - d{t) - [{ey- l)^y\>^}g{t,y)dy. 
Jr 

4.1. Time-inhomogeneous chains. To approximate the barrier option prices driven by a given 
time-inhomogeneous Markov process we again start by constructing an appropriate Markov chain. 
We build an approximating time-inhomogeneous chain with a generator that is piecewise constant 
in time. Given a partition T = {7i}"^Q with Tq = < Ti < ... < T„ = T of [0, T], we assume that 
the chain X on the state-space G evolves according to the generator A*^*) during the time-interval 



proof of these facts is given in Appendix 1X1 
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[Tj_i,Tj) where the matrices A^*^ are chosen so as to approximate well the infinitesimal generators 
Ct^ . Thus X has a time-dependent generator given by 



(4.7) Ai:=^A«I[^^_,,^^)(t), t>0. 

i=l 

Also we take the short rate r*^") to be piecewise constant 

n 

(4.8) r^''\t) = Y,nl[T,_,,Ti){t), for t >0 and >0, i = l,...,n. 

i=l 

As a straightforward consequence of Theorem [T] we have the following result for the first-passage 
probabilities of the chain X. 

Corollary 1. For r^") given as above and any functions : G ^ M and ^ : G — > R with ip{x) = 
for X G E, the following equalities hold 

(4.9) B,[cl>{XT)I{r>T}] = (exp(AriAW)---exp(Ar„A("))(/>) (x), X G G, 

(4.10) E^[e-fo''"Ht)dt^^Xr)I{r<T}] = (exp(AriA«)...exp(Ar„A(::))V^)(x), x G G, 
where ATj := Tj — Tj_i, and A*^') and Ar*^ are defined as in ()3.5p and ()3.4p . 

5. Construction of the Markov chain 

In this section we construct the approximating Markov chains used in the algorithm introduced 
in this paper. In Section [5.11 we review the discretization of diffusion processes and in Sections 15.21 
and 15.31 we consider the Markov processes with discontinuous paths. In Section [5.41 we describe the 
algorithm for time-inhomogeneous Markov processes with and without jumps. 

5.1. Diffusions. Let S = {St)t>o be an asset price process which evolves under a risk-neutral 
measure according to the SDE in ()2.7p of ExamplelH For a given finite state-space E of the diffusion 
S (and a given barrier contract) we construct a non-uniform state-space G with elements using 
the algorithm in Appendix [Bl Define the sets 

dG:={xi,XN} and G° := G\aG, 

where the "boundary" dG consist of the smallest (i.e. xi) and largest (i.e. xn) elements in G and 
the "interior" G° is the complement of the boundary. 

The next task is to construct the approximating continuous-time Markov chain X = {Xt)t>o 
on the state-space G by specifying its generator matrix A in such a way that the first and the 
second instantaneous moments of the processes 5* and X coincide on the set G°. In other words 
the following condition needs to be satisfied: 

(5.1) [{SAt - SoY] = Exo [{XAt - XoY] + o{At), for Xq e G° j e {1, 2}, 

where = Xq. Well-known facts in the theory of diffusion processes and continuous-time Markov 
chains imply that the coefficients of the generator matrix A must satisfy the following system for 
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each X £ G°: 



(5.2) 



and A{x,y)>0 Vy G G\{x} 



(5.3) 



^A{x,y){y - x) 



IX, 



(5.4) 




For x S 5G we impose an absorbing boundary condition: A(x, y) = for all y G G. The condition 
in (j5.2p ensures that A is a generator matrix and equation (j5.3p implies that the chain X drifts 
locally at the same rate as the diffusion S. 

In applications it is typically possible to find a tri-diagonal generator matrix A that satisfies 
the linear system in (15. 2p - (15. 4p . In terms of the process X this implies that at any moment 
in time the chain can only jump to neighbouring states. By setting the top and bottom rows of 
the matrix A to zero, we make the states in the set 9G absorbing. This behaviour of the chain 
X clearly differs from the dynamics of the diffusion S in the neighbourhood of the set (9G. We 
therefore have to choose the boundary states far enough that the laws of the processes X and S 
are close to each other during the finite time interval of interest. Such a choice is computationally 
feasible because the state-space G is non-uniform. In practical applications we can ensure easily 
that the accumulation of probability mass in the states of 5G is negligible during the time interval 
of interest. 

Note that the Markov chain approximation X of the diffusion S defined by the linear system 
in (j5.2p - (j5.4p is by no means the only viable alternative. One could in principle produce more 
accurate results by matching higher instantaneous moments of the two processes. However for the 
sake of simplicity we do not pursue this idea further at this stage. 

In Sections 17. II and 17. 21 we the compare numerical results of our algorithm with the corresponding 
results for the geometric Brownian motion (reported in [21] and [31] ) and the CEV process (reported 
in [H]) respectively. 

5.2. Levy subordinated diffusions. A common way of building models with jumps is by time- 
changing diffusions with a Levy subordinator (see e.g. Sections 17.31 and l7.4p . If the jump process 
has this special structure, its generator is closely related to the generator of the diffusion. In this 
section we recall how this relationship can be exploited to obtain a continuous-time Markov chain 
approximation of the jump process that admits this special structure. 

Let S' be a time-homogenous diffusion governed by SDE (|2.7p that satisfies the conditions in 
Example and let Z he a Levy subordinator (i.e. a Levy process with non-decreasing paths that 
starts at 0) independent of S'. In this section we assume that the asset price process is modelled 
by the process S = {St}t>o, where 



It follows from the Phillips theorem (see [481 Thm. 32.1] or [44J) that S" is a Markov process that 
satisfies the Feller property. It follows by conditioning on the independent subordinator Z that the 



(5.5) 



e^*5', 



for some /i G M. 
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discounted asset price process S = {e'^^^'^^^ St}t>o is a martingale if and only if [St] = e^^^'^^^So 
for alH > and So G E, which is equivalent to 

(5.6) fi + ^zi-l) = r - d 

where V'z is the Laplace exponent of the subordinator Z (i.e. Eo[e~"^*] = e*'^^^"^), 7 is the drift in 
SDE (j2.7p and r, d are the interest rate and the dividend yield respectively. The identity in (j5.6p 
follows by conditioning on the independent subordinator Z and using the fact E5(, [Sf] = Soe^^ for 
all t > and Sq £ E. If (|5.6p is satisfied, we can use as a model for the risky asset under a 
risk-neutral measure without introducing arbitrage into the market. 

We will now use the structure of the process S given by subordination ()5.5p to define a generator 
matrix for the Markov chain that approximates the asset price process S. Recall first that the 
generator Q' of the diffusion S' can be expressed as 

{g'f){x):=jx{Vf){x) + ^x^a{xf{Af)ix), for / G ^^(E), 

where (A/)(x) = f"{x) and {W f){x) = f'{x). Recall that by the Phillips theorem (see [48', Thm. 
32.1] or [H]) the infinitesimal generator Q'^ of the process {S'^^}t>o can be expressed as 

(5.7) g'^ = ^z{-G'), 

where is the Laplace exponent of the Levy subordinator Z. Formula (j5.7p is a formal identity 
but can be interpreted in terms of functional calculus (see the references in the paragraph following 
Lemma [T] in Section [2]) . Representation (|5.7p and the general theory of Markov processes imply 
that the infinitesimal generator Q of the process S is given by the following expression 

(5.8) (g/)(x) = MV/)(x) + (g^/)(x) for any f € C^iK). 

In order to construct the approximating Markov chain we fix a finite state-space G C E (using 
the algorithm described in Appendix |Bj) and define a tri-diagonal generator matrix A' by the 
linear system in (I5.2p ~ (l5.4p . the right-hand side of which is given by the drift coefficient 7 and the 
diffusion coefficient a of SDE (12. 7p satisfied by the process S' . It is clear that the Markov chain 
X' , which corresponds to the generator A', is by construction a Markov process that approximates 
the diffusion S' in the sense of Section 15.11 

We now define a generator matrix A^ by applying the Laplace exponent ipz of the Levy subordi- 
nator to the matrix —A'. Since the semigroup of the chain X' is generated by a bounded operator 
A', Phillips' theorem (see [IHl Thm. 32.1] or [Hj) implies that the matrix A'^ := tpzi—-^') is a gen- 
erator matrix of the (time-changed) Markov chain {X'^^}t>o on the state-space G. It is therefore 
natural to take the chain generated by A^ as an approximation for the process {S'^^}t>o- Note 
that in order to compute the generator A^ we first obtain the decomposition A' = UDU~^, where 
D is a diagonal matrix with eigenvalues equal to those of A' and U is an invertable matrix, and 
define A^ := U4'z{—D)U^^ , where ipzi—D) is a diagonal matrix with eigenvalues tpzi—^) when 
A ranges over the spectrum of A'. In practice the matrix A' always has a diagonal decomposition 
because the set of all square matrices that do not possess it is of codimension one in the space of 
all square matrices and therefore has Lebesgue measure zero. 
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In order to obtain a matrix that approximates the operator given in ()5.8p . we have to add drift 
with state-dependent intensity to the Markov chain generated by the matrix A^. We define a 
tri-diagonal generator matrix A^ that satisfies the hnear system 

(5.9) J];A^(x,y)(y-x) = {r - d)x -Y,^z{^,z){z - x) Vx e 

and A^(x,y) = for all x S 9G, y E G. The generator matrix A that is used to approximate the 
operator in ()5.8p can now be defined as 

A := A'^ + A/,. 

Intuitively we add to A^ a matrix A^ which has in each row at most one non-zero element of the 
diagonal. In other words, if the right-hand side of (|5.9p is positive (resp. negative) we increase the 
value of the element on the superdiagonal (resp. subdiagonal) in the corresponding row of A'^, thus 
increasing the overall intensity of the approximating chain to jump in the required direction. This 
construction ensures that the chain X drifts locally at the same rate as the original process S. As 
in the diffusion case, the top and bottom rows of the generator A are set to be equal to zero. 

5.3. Jump processes with state-dependent characteristics. In this section we consider Markov 
processes with state-dependent characteristics. The basic idea is to define the approximating 
Markov chain so that the first two instantaneous local moments of the original process are matched, 
i.e. so that condition ()5.ip holds. We now describe the procedure in more detail. 

The asset price process S is in this section assumed to be a non-negative Markov process with 
the generator in (j2.10p with a possibly state-dependent volatility cj{x) and jump measure v{x^dy). 
The construction of the generator matrix A of the approximating Markov chain X is now carried 
out in two steps: we first define the jump matrix Aj, which corresponds to the discretization of 
the jump measure and then characterize a tri-diagonal generator matrix Ac by stipulating that 
the Markov chain X with the generator A = Aj -|- Ac has the same instantaneous local moments 
as the process S. 

We start by building the state-space G C E with N elements using the algorithm in Appendix iBl 
For any given x E G° we transform the set G into the set G^ C (— 1, oo) defined by 

G^ := 1^ - 1 : z E g} . 

It is clear that the set G^ consists of the relative jump sizes of the approximating chain X. It 
is therefore natural to define the x-th row of the jump part Aj of the generator of X using the 
discretization of the jump measure z^(x, dy) on the set G^,-. In particular let {yj : i = 1, . . . , A^} = G^;, 
where yi < y^+i for alH = 1, . . . , — 1, yo ■= —1 and yAr+i := oo and define a function 



(5.10) 



Ox : GxU {yo} ^ [-l,oo] by a^iyo) = -I, ax{yN) = oo, 
axiVi) ^ {yi,yi+i) for i E {1, . . . , - 1}. 
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A possible natural choice for ax{yi) would be the middle point of the interval We can 

now define the jump part of the generator as 

(5.11) Kj{x,x{l + yi)) := v {x,{ax{yi-i),ax{yi))) where z G {1, . . . , iV} end / 0, 

(5.12) Aj(x,x) = - ^^(^'^)- 

Note that the formula in (j5.11|) simply allocates the jump intensity of a relative jump of the 
process S in the interval {oLx{yi-i) , axiVi)) to the jump intensity of a relative jump of size yi of 
the approximating chain X. Since ^{x^ dy) is a measure, the expression in (|5.12p ensures that the 
matrix Aj is a generator matrix. For x € dG we set Aj{x,y) := for all y £ G. This completes 
the first step in the construction of the generator matrix A of the chain X. 

In the second step we match the first and second instantaneous moments of the asset price process 
S. In other words the chain X must satisfy condition (j5.ip for all starting states Xq = Sq £ G°. 
Note that condition (|5.ip implicitly assumes that the second instantaneous moment of S exists. 
This is the case if the jump measure satisfies the following condition 

/oo 
y'^u{x,dy) < oo Vx G E 

which we now assume to hold. 

The task now is to find a tri-diagonal generator matrix Ac such that the chain generated by 
the sum Ac + Aj satisfes (|5.ip . The tri-diagonal matrix Ac therefore has to satisfy the following 
conditions 

(5.14) ^Ac(3;,z) = and Ac(x,z)>0 Vz € G\{x}, 

(5.15) 'YAc{x,z){z-x) = {r - d)x - 'Y Aj{x,z'){z' - x), 

zeG z'eG 



(5.16) > Ac{x,z){z-xf = X 



zGG 



a (xf 



/oo "I 
y'^u{x,dy) - ^ Aj{x,z'){z' 
"1 J '.tr-(n 



X? 



'eG 

for all X G G°, where r, d are the instantaneous interest rate and dividend yield respectively and a is 
the local volatility function in (|2.1Up . The right-hand side of equation (|5.15p is the difference of the 
risk-neutral drift and the drift induced by the presence of jumps. Similarly the right-hand side of 
the linear equation in (I5.16P consists of the difference of the instantaneous second moments of the 
asset price process S (computed directly from its generator ()2.10p ) and the chain that corresponds 
to the jump generator Aj. As usual we assume the absorbing boundary condition Ac(x, y) = for 
all X G dG, y G G. 

The linear system in (j5.14p - (|5.16p can typically be satisfied by a tri-diagonal generator matrix 
Ac, because for a given jump measure i^(x,dy) we can choose the function in ()5.10p so that 
the right-hand side of (I5.16P is positive. Once we find Ac we define the generator matrix of the 
approximating chain X by 

A := Ac + Aj. 
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Note that condition ()5.16p implies that the infinitesimal drift of the chain X is the same as that 
of the Markov process S. By stipulating in ()5.16p that the instantaneous variances of S and X 
coincide, we improve the quality of the approximation and hence provide more accurate prices of 
the contingent claims. Finally we note that the absorbing boundary condition on the generator 
matrix A does not influence the law of the approximating process if the boundary points are chosen 
to be far enough from the current spot price. 

Note that the algorithm outlined above can be applied to any Feller process with state-dependent 
characteristics and finite second instantaneous moments, including the processes with infinite ac- 
tivity and infinite variation. We will illustrate this algorithm in Section 17.51 in the case of a local 
Levy process. 

5.4. Time-inhomogeneous Markov processes. In the case of diffusion (|4.2|) with time-dependent 
coefficients S and 7, we have to find an approximate generator for each time step and then apply 
the formula from Corollary [T] to obtain the first-passage probabilities and barrier option prices. The 
algorithm is in this case a straightforward generalisation of the procedure described in Section 15.11 
The numerical results are presented in Section I7.6[ 

Assume now that the asset price process S is an exponential additive process given by (j4.3p 
with generator (j4.6p that depends on time (i.e. the process S has both jumps and time-dependent 
characteristics). It is clear that during a short time interval of constancy a modification of the 
algorithm described in Section 15.31 can be applied to obtain the generator of the approximating 
chain. Once the algorithm has been applied for each of the short time intervals. Corollary [T] can be 
used to find the first-passage probabilities and the corresponding barrier option prices. 



6. Weak convergence 

We next turn to the question of how to construct a sequence of finite-state continuous-time 
Markov chains X^") such that the corresponding expected payoffs approximate barrier option prices 
under a given Feller price process S = {St}t>o- For X^") to replicate as closely as possible the 
dynamics of S one chooses the generator matrix A*^"-* with the corresponding state-space G^"^ such 
that it is uniformly close to the infinitesimal generator C of S", in the sense that the distance e„(/) 
between the generators is small for a sufficiently large class P of regular test functions /, where 

€„(/) := max AW/„(:e) - Cf{x) 

and /„ = /|(g(n) is the restriction of / to G^"). More specifically, if e„(/) tends to zero as n tends 
to infinity for / in the class P, then the sequence of processes (X("))„gpsi converges weakly to 
the process S. This weak convergence on the level of the process implies in particular that the 
marginal distributions of X^"^ will converge to those of S, and therefore the values of European 
options converge, that is, 

E,[/(4"))] ^ E,[/(St)] 
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for E E, maturity T > 0, and continuous bounded functions /. In practice the boundedness of 
the payoff / is not restrictive as it is always possible to consider the truncation / A M for constants 
M large enough without losing noticeable accuracy. 

The payoff of the barrier option can be described in terms of the first-passage time of S and 
the position of S at that moment, which are both functionals of the path {St}t>o- For the weak 
convergence of X^") to S to carry over to convergence of barrier-type payoffs, continuity (in the 
Skorokhod topology) is required of these two functionals, which is guaranteed to hold under As- 
sumption [2j In view of the fact that the payoff of a barrier option is typically a discontinuous 
function, an additional condition is needed to ensure the convergence of the barrier option prices; 
we will assume that i and u are such that 



(6.1) 



P,iSTG{i,u})=0. 



Most models used in mathematical finance satisfy this condition. Even if (16. ip is not satisfied, this 
does not constitute a limitation in practice, since for any given process S the condition is satisfied 
for all but countably many pairs {i,u). The statement of the convergence is made precise in the 
following theorem, the proof of which is in Appendix [Xl 

Theorem 2. Let S be a Feller process with state-space E and infinitesimal generator C satisfying 
Assumptionl^ and X^"'^ a sequence of Markov chains with generator matrices A^") such that 



as n 



oo 



(6.2) €n{f) ^ 

for any function f in a corj^ of C If (16. ip holds, then, as n 



oo. 



E, 



9 { 4"^ 



E, 



^{r(")>T} 



for any bounded continuous functions g,h 



E, [g{ST)I{r^>T}] , 

E, [e— ^/i(5.JI{,^<r}] , 



where r 



(n) 



inf{t > : ^ A}. 



For the time-inhomogeneous case an analogous convergence result holds true. To approximate a 
given time-inhomogeneous Markov process {St}t£[o,T]^ the sequence of time-inhomogeneous gener- 
ator matrices A^") of the form (j4.7p . defined on the time-space grids T*^") x G^"^ needs to be close 
to the (space-time) generator C of S (see Section H] for the definition of £'), where the distance is 
measured by 



enif) := niax 



AI""^ ft,n(.x) - Ctftix) 



for / in a core of C. Here ft denotes the function ft:x>-^ fit^x) and ft^n = ftl^i") restriction 
of ft to G("). Further, in the case of the rebate options a given short rate r{t) also has to be 
approximated by an appropriately chosen piecewise constant function r^") as in (|4.8p . 



A core is a dense subset of Co(E) such that the set {(A — £)/ : / G Co(E)} is dense in Co(E) for some A > 0. 
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Corollary 2. Let {D,Y) be the Feller process defined in ()4.1|) with state-space E' = [0,T] x E and 
infinitesimal generator C , and let X^^^ be a sequence of time-inhomogeneous Markov chains with 
generator matrices (A^"\t S T^")) such that 

(6.3) ^0 as n ^ oo, for any function f in a core of CJ . 

If r{t) is continuous and S satisfies (16. ip and AssumptionlM (with replaced by Pq^x), then, as 



n 



oo, 



5(4"^ 



{r("'>T} 



E, [<7(5t)I{.^>t}] 



E, 



,(n) j I{r(")<T} 



E, 



TA)HrA<T] 



for any bounded continuous functions g,h -.E, 



where stopping times are as in Theorem\^ 



7. Numerical examples 

In this section we are going to examine numerically the behaviour of our algorithm in a variety 
of contexts. 

7.1. Geometric Brownian motion. The model is given by SDE (|2.7p where the volatility func- 
tion ct{x) = (To is constant and the drift equals ^ = r — d, where r is the risk free rate and d 
is the dividend yield. We now compare our algorithm (MG), based on the Markov generator of 
the approximating chain X, with the results obtained in Geman and Yor [2T| and Kunitomo and 
Ikeda |34| . The numerical results are contained in Table [TJ 



ctq = 0.2, r = 0.02 
K = 2, £ = 1.5, u = 2.5 


erg = 0.5, r = 0.05 
K = 2, £ = 1.5, w = 3 


CTo = 0.5, r = 0.05 
K = 1.75, £ = 1, u = 3 


GY 


KI 


MG 


GY 


KI 


MG 


GY 


KI 


MG 


0.0411 


0.041089 


0.041082 


0.0178 


0.017856 


0.017856 


0.07615 


0.076172 


0.076165 



Table 1. The comparison of double barrier option prices obtained in and 01] in the case of 
geometric Brownian motion. The model is given by SDE (|2.7p with the constant volatility function 
cr{x) := (Jo and drift 'y — r ~ d, where the interest rate r is given in the table and the dividend 
yield equals d = 0. The asset price process S starts at So = 2 and the maturity in all the cases is 
T = 1 year. The state-space of the approximating chain is defined by the algorithm in Appendix IbI 
and the parameters TV = 200, minS = 0.2,maxS = 10, = 100, 5^ = l,gf = 10, fff = 10, = 
l,g^ = 100. The computation for the pricing of the barrier products takes about 0.03 seconds (on 
Intel(R) Xeon(R) CPU E5430 @ 2.66GHz) for each of the parameter choices in this table. 



7.2. CEV process. The model is given by SDE (|2.7p where the volatility function takes the form 
a{x) := ao{x/So)^ {Sq is the starting value of the process) and the drift equals 'y = r — d where 
r is the risk free rate and d denotes the dividend yield. Table [2] contains the numerical results of 
our algorithm (MG) and compares them to the results obtained in Davydov and Linetsky [Hj. See 
also [l3] for the implementation of the pricing algorithm in Matlab that produced these results. 
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CEV 


13 = 


(3 = -0.5 


13 = -I 


u ( K 


MG: n = 2 


MG: n = 3 


DL 


MG: n = 2 


MG: n = 3 


DL 


MG: n = 2 


MG: n = 3 


DL 


120 90 95 
120 90 100 
120 90 105 


1.7038 
0.9703 
0.4417 


1.7038 
0.9703 
0.4418 


1.7039 
0.9703 
0.4418 


1.8805 
1.0957 
0.5124 


1.8805 
1.09-57 
0.5125 


1.8805 
1.0958 
0.5126 


2.0799 
1.2382 
0.5943 


2.0800 
1.2383 
0.5944 


2.0800 
1.2383 
0.5945 




/3 = -2 


/3 = -3 


/3 = -4 


u I K 


MG: n = 2 


MG: n = 3 


DL 


MG: n = 2 


MG: n = 3 


DL 


MG: n = 2 


MG: n = 3 


DL 


120 90 95 
120 90 100 
120 90 105 


2.5527 
1.5797 
0.7958 


2.5528 
1.5798 
0.7960 


2.5529 
1.5799 
0.7960 


3.1292 
2.0019 
1.0532 


3.1294 
2.0021 
1.0534 


3.1295 
2.0022 
1.0535 


3.8084 
2.5055 
1.3693 


3.8087 
2.5058 
1.3696 


3.8088 
2.5059 
1.3696 



Table 2. The comparison of double barrier prices obtained in [TJ using the spectral decompo- 
sition of the CEV process and our algorithm based on the Markov generator. The model is given 
by SDE (|2.7p with the volatility function a{x) := ao{x/So)^ and drift j = r — d, where r is the 
interest rate and d is the dividend yield. The model parameters are as follows: So = 100, ao = 0.25, 
r = 0.1, d = and all options in the table expire at maturity T = 0.5. The state-space of the 
approximating chain is defined by the algorithm in Appendix [B] and the parameters N = 2" ■ 100 
(where n = 2,3), minS = 0.1,maxS = 190, gt = 100, = l,gf = 10, = 10, = l,g^ = 100. 
The computation time for the pricing of all the barrier options takes about one second in the case 
n = 2 and about ten seconds when n = 3 on the same hardware as in Table [T] for each of the six 
models in this table. 



7.3. GH process. The price process S in this example is assumed to be an exponential Levy 
process given by (12. 8p where L is a generalised hyperbolic (GH) Levy process with distribution at 
time 1 given by the characteristic function 

(7.1) «l>z.,(n) := Eo e-^^ = e^M , .o^. ,2 ) ^ , X — 



where m, A G M, a, 5 > 0, < < a. 



The function K\ is known as the modified Bessel function of the second kind (see p!] for the precise 
definition of K\). The corresponding distribution is called generalised hyperbolic and is denoted 
by GH(A, a, /3, (5, m). The class of distributions described by the characteristic function in ()7.ip 
was introduced into the mathematical finance literature by Eberlein and Keller [18] and has been 
studied extensively ever since (see e.g. Prause [35] and the references therein). Since L has to 
satisfy the exponential moment condition (|2.9p . we have to restrict the parameter space in the 
following way (see for example Lemma 1.13 in [45j for details): 



(7.2) m, A e M, a,5 > and 



1 

+ -<«. 



Furthermore it is clear from ()2.8p that the value of the parameter m in (17. ip has no bearing on the 
model S and can without loss of generality be taken equal to zero. 
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Generalised hyperbolic Levy process 
a = 26.4, [i = -0.53, 5 = 0.0034, A = -0.5 


n = 1 


n = 2 


n = 3 


n = 4 


n = 5 


KO Call: K = 100, ^ = 97, m = 105 


1.0193 


1.0190 


1.0188 


1.0187 


1.0187 


Double-no-touch: I = 97, u = 105 


0.9704 


0.9711 


0.9714 


0.9716 


0.9716 



Table 3. The prices of the double barrier knock-out call option and the double-no-touch option 
in the generalised hyperbolic Levy model. The parameter values for the Levy process L axe taken 
from [45], page 64, Table 2.27, and are given in the table. The strike in the case of the call option 
is given hy K = 100 and the barrier levels for both derivatives are equal to ^ = 97, u = 105. 
We take = 100, r = 0.01, d = and T = 1. The risk-neutral drift ^ in (f^ takes the 
form /i = r — d — tlizi—P — 1/2) ~ —0.0099. The state-space of the approximating chain is 
defined by the algorithm in Appendix [Bl and the parameters TV = 2" • 100 (where n = 2, . . . , 5), 
minS = 50,maxS = 150, gt = 100, = l,gf = 10, = 10, = l,gu = 100. The approximate 
computation time of the pricing algorithm (on the same hardware as in Table [1} for the case n — 3 
is 19 seconds. 



Let Z he a generalised inverse Gaussian (GIG) subordinator with the Laplace exponent tpz given 

by 



/ -s\/2Kx( Jb{a + 2u) 
(7.3) Eofe-^^i] = e'^^(")= — ^ where A G M, a, 6 > 0. 



The distribution of Zi is called generalised inverse Gaussian and denoted by GIG(A, a, b). It follows 
from expressions ()7.ip and (j7.3p that if we choose Zi ~ GIG(A,q^ — (3'^,S'^) where the parameters 
A, a, 13, 6 satisfy condition (|7.2p . then the equality 

Eo[ei"^i] = exp(V^z(nV2-i/3n)) 

holds for all n G M. Therefore L has the same law as the process {(3Zt + Wzt}t>o where W is 
a Brownian motion independent of Z. Hence the asset price process S obviously posseses the 
structure in (jS.Sp where the diffusion S' is a geometric Brownian motion and the drift equals 
^ = r — d — V'z(— /? — 1/2)- We can therefore apply the algorithm described in Section [5^21 Table[3] 
contains the numerical results. 

7.4. The CGMY/KoBoL process. In this section we assume that the price S is again an ex- 
ponential Levy process given by ()2.8p . where the Levy density of the process L is given by the 
formula 

/ e-G\y\ Q-My \ 

(7.4) k{y):=C\ I|j,<o} + I{y>o} 1 , where M>1,G>0, C>0,y<2. 

The inequality y < 2 is induced by the integrability condition on the Levy measure at zero and 
the condition M > 1 implies the exponential moment condition (j2.9p . 

Madan and Yor [52] show that there exists a Levy subordinator Z with the Laplace exponent 
V'z given by 

(7.5) Eo[e-"^*] = e*^^(") = exp {tCT{-Y)[2r{uf cos(r/(n)y) - - G^]) , u > 
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Spot 


Knock-out put 


Double-no touch 




Spot 


Knock-out put 


Double-no touch 


% 


BL 


MC 


MG 


BL 


MC 


MG 


% 


BL 


MC 


MG 


BL 


MC 


MG 


82 


302.28 


301.88 


301.07 


0.5778 


0.5768 


0.5757 




101 


64.07 


64.56 


64.53 


0.9481 


0.9487 


0.9483 


85 


370.38 


370.98 


370.38 


0.8009 


0.8004 


0.8004 




104 


36.96 


37.13 


37.18 


0.9351 


0.9357 


0.9352 


88 


341.35 


341.73 


341.78 


0.8881 


0.8887 


0.8880 




107 


22.73 


22.77 


22.84 


0.9112 


0.9115 


0.9113 


91 


279.86 


280.71 


280.41 


0.9280 


0.9288 


0.9280 




110 


14.60 


14.58 


14.65 


0.8708 


0.8703 


0.8709 


94 


207.71 


208.48 


208.30 


0.9464 


0.9463 


0.9465 




113 


9.61 


9.58 


9.64 


0.8020 


0.8026 


0.8019 


97 


136.63 


137.23 


137.24 


0.9527 


0.9531 


0.9529 




116 


6.30 


6.27 


6.32 


0.6771 


0.6773 


0.6767 


100 


78.19 


78.67 


78.74 


0.9506 


0.9508 


0.9507 




119 


3.52 


3.54 


3.54 


0.4049 


0.4063 


0.4049 



Table 4. Barrier option prices under the CGMY model. The first column contains the spot 
price as percentage of 3500. The CGMY parameters are C = 1, G = 9, M = 8, F = 0.5. The 
resulting risk-neutral drift is /i = r — d — tpzi—0 — 1/2) ~ —0.0423. Option parameters K = 3500 
(strike of the put), I = 2800, u = 4200, r = 0.03, d = 0, T = 0.1. The columns BL and MG report 
the results obtained by Boyarchenko and Levendorskii [5] and by the Markov generator algorithm 
(with A*' = 800 points) respectively. The column MC gives the results of the Monte Carlo pricing 
algorithm given in ^5. . It takes about 22 seconds to run the MG algorithm for each starting spot 
price on the same hardware as in Table [1] 

where T denotes the Gamma function and the functions r, 7] are given by the formulae 
r{u) := V2u + GM and r/(n) := arctan 

Note that for u £ [-GM/2, 9^/2) the function 7] in formula (17. 5|) takes purely imaginary values 
which are mapped by the cosine function into the real numbers. Furthermore it is shown in [32] that 
the Levy process L, given by the Levy density ()7.4p . has the same law as the process {Wzt+OZt}t>Q 
where 9 := {G—M)/2 and the Brownian motion W is independent of the subordinator Z. Therefore 
we have the identity Eo[exp(Li)] = exp(V'z(— 6*— 1/2)) which (together with definition (j2.8p ) implies 
that the process S has the same law as the process given by (|5.5|) where S' is a geometric Brownian 
motion and the drift ^ satisfies ^ = r — d — tpzi—G — 1/2) by (|5.6p . We can therefore apply 
the algorithm described in Section [5.21 We compare our numerical results with those obtained in 
Boyarchenko and Levendorskii (see Table H]). 

7.5. Local Levy model. A Markov process S with state-dependent characteristics, which starts 
at So G E and has double-exponential jumps can be specified by the following Levy driven SDE 

{r-d- XC{St-/Sof)dt + {St-/SofdLt, where 

Nt 

(ToWt + (e^» - l) , do G (0, oo) and /? G M. 

1=1 

The special case of this model for /? = was introduced into the mathematical finance literature by 
Kou [32]. The random variables ifj, i G N, are independent of both the Brownian motion W and 
the Poisson process N with intensity A > and are distributed according to the double exponential 



I G + M 



(7.6) 



(7.7) 



dSt^ 

St- 



Lt 
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density 

(7.8) fK{k) = priie-'''%k>o} + {l-p)me'^'%k<o}, where 

771 > 1, ?72 > and p G [0, 1]. 

The parameter ^ is given by 

C:=Ere^^-ll=^ + ii^-l. 

r?i - 1 ??2 + 1 

It is clear that the model described by (j7.6p and (j7.7p has a generator of the form given in ()2.10|) 
with a{x) = (Tq{x/SqY and jjL{x) = \Q{x / Sq)^ . The jump measure v{x,dy) in representation (|2.1U|) 
is supported in (—1, 00) and in our case by ()7.8p takes the explicit form 

(7.9) u{xAy) = (x/5o)^A[pr?i(y + l)-i-''iI|,>o} + (l-p)^2(y + l)''^-'l{-i<j,<o}]d?/. 

Note that the drift ^(x) and the jump measure v{x^dy) satisfy the condition in (j2.1ip . 

Representation (j7.9p of the jump measure z/(x, dy) of the asset price process S can now be used 
to construct the "jump" generator Aj defined in equations (jS.lip and (j5.12p . Furthermore it is 
clear from (17. 9p that the instantaneous variance term caused by the jumps of the process S in (I5.16P 
is of the form 

yMxAy) = {x/Sof2\(- +^lZP_\ if >2. 

-1 V(m - l)(f?l - 2) (7?2 + l)(?/2+2)y 

Note that condition (I5.13P in the context of the present model is satisfied if ryi > 2 which is typically 
true in applications. The linear system in (|5.14p ~ (j5.16p can now be solved. The numerical results 
of the MG algorithm applied to an up-and-in call option are contained in Table O where they are 
compared (in the case /3 = 0) with the corresponding results of Kou and Wang [HS] (see Table 3 
in EHI). 



7.6. Time dependent jump-diffusion. Let L be a compound Poisson process with intensity A 
and jumps of the form (e^ — 1)) where K is a. normal random variable with mean m and variance 
s^. Consider an asset price process S which starts at S'o G IE and is given by the SDE 

A Q 

(7.10) = {r{t)-K)dt + ^{t,St.)mt + dLt, where 
Jt- 

C := E [e^ - 1] = e'"+^'/2 _ ^ 

and the Brownian motion W is independent of L. The volatility function S(t, x) and the instanta- 
neous interest rate r{t) are given by 

(7.11) S(t,x) := v{t){x/Sof, where v{t) := 9 + {a^ - 9)e-''\ 

(7.12) r{t) := tq + rie''°*. 

The constant v{0) = ctq is the starting value of the volatility in our model while the constant 9 
represents the level to which the function v tends in time. The process S can be viewed as a local 
volatility Merton jump diffusion. Such processes were used in [3] and [25] for calibration purposes. 
The parametric form of the time-dependent local volatility function has been used in practice and 
was taken from HO], Section 10.7. 
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Local Levy model 


TV /T f '^ ~\T A r\r\ 

MG: N = 400 


MG: N = 800 


MG: N = 1200 


KW 


/3 = 


A = 3 


10.0528 


10.0530 


10.0530 


10.05307 




A = 0.01 


9.2768 


9.2771 


9.2772 


9.27724 


/?=-! 


A = 3 


9.7685 


9.7688 


9.7688 


N/A 




A = 0.01 


8.9572 


8.9575 


8.9575 


N/A 


/3 = -3 


A = 3 


9.0185 


9.0187 


9.0188 


N/A 




A = 0.01 


8.0855 


8.0858 


8.0858 


N/A 



Table 5. The model 5" is defined in (|7.6|l and (|7.7p . The parameters are given by 5*0 = 100, 
r = 5%, d = 0, cro = 0.2, p = 0.3, l/r?i = 0.02, l/r?2 = 0.04, /3 G {0, -1, -3} and A £ {0.01, 3}. 
The strike in the up-and-in call option \& K = 100 and the upper barrier u = 120 while time to 
maturity T — 1. The column KW denotes the results of Kou and Wang (see [33], Table 3) while 
MG denotes the algorithm based on the Markov generator. The state-space of the approximating 
chain is defined by the algorithm analogous to the one in AppendixjB] adapted in an obvious way to 
a single barrier contract. Its size is = n ■ 400 for n = 1, 2, 3. The computation time is about half 
a second for n = 1, five seconds for n — 2 and seventeen seconds for n = 3 on the same hardware 
as in Table [1] We also ran the algorithm for A'' = 1600 and obtained identical results (up to four 
decimals) as the ones in the column — 1200. 



The Levy measure of L is of the form 

^7 1Q^ fA \ ^ f {log{l + y) -m)^ \ dy 

(7.13) u{dy) = ^Tl^exp^ J yG(-l,oo). 

If $ denotes the standard normal cdf we find that if — 1 < a < 6 then 

f i^{dy) = A[$((log(l + 6) -m)/s) -$((log(l + a) -m)/s)] and 

J a 

/oo 
y2^(dy) = A [e2(-+«') - 2e"^+^'/2 ^ 1 . 

We can therefore apply directly the discretization algorithm described in Section [5.31 to obtain the 
generator of the approximating chain X during each time interval of constancy. An easy application 
of Corollary [T] yields the results in Table El 

7.7. VG-Sato process. Here we implement the algorithm described in Section f5. 41 in the case of 
the exponential Sato process, which was introduced into financial modelling by Carr et al. [lOj . 
Recall that a Sato process is an additive process X = {Xt)t>o, which is self-similar (i.e. Xt ~ f^Xi 
for some constant 7 > and all t > 0) and whose law at time one is self decomposable. In Theorem 1 
of [TU] it is proved that the characteristic function of Xt is of the form 

(7.14) $x(n,t) = Eo[e^^^*] = exp (^^ (e^"^' - l) ^^^^dy) 
where the function /i : R — > satisfies the conditions 

f h(x) 

h(±xi) > h(±X2) for all < xi < 2:2 and / min{x^, l}-j — r-dx < 00. 

Jr m 
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Time dependant jump-diffusion 


^ in 


/i — iUU 


N = 400 


Rebate: i = 90,u = 120 
Double-no-touch: £ = 90, it = 120 


0.4577 

U. loDo 


0.4574 

U. 10 / z 


N = 800 


Rebate: £ = 90,u = 120 
Double-no-touch: £ = 90,u = 120 


0.4570 

n 1 '^fiv 
U.loD / 


0.4567 


N = 1200 


Rebate: e = 90,u = 120 
Double-no-touch: £ = 90,u = 120 


0.4568 

U.ioo ( 


0.4565 

n 1 QVo 
U.io ( / 


N = 1600 


Rebate: £ = 90, m = 120 
Double-no-touch: £ = 90,u = 120 


0.4568 
0.1367 


0.4565 
0.1372 



Table 6. The model is jump-diffusion (|7.10p with volatility function (|7.11|) . time-dependent 
drift (|7.12p and normally distributed jumps (I7.13p . The constants in the volatility and drift function 
are given by So = 100, ao = 0.25,6 = 0.2, fc = 10,(3 = -2 and ro = 0.01, ri = 0.09, ao = -1. The 
function r{t) in (I7.12p represents the instantaneous time-dependent interest rate (the dividend yield 
equals d{t) = 0). The jump parameters A — 0.1008, m = —0.9144 and s = 0.4367 are taken from 
the calibrated model in [53]. The state-space of the approximating chain is defined by the algorithm 
in Appendix [B] with the same parameters as in Tabled A'^ is the number of states and n stands for 
the number of time steps. The table contains the prices of the double-no-touch option that pays 
one at expiry if the barriers £ and u have not been breached, and a rebate option (i.e. American 
knock-in option) that pays one when the spot leaves the corridor between the barriers. All options 
expire in T = 0.5 years. The algorithm was also run for n = 500 time steps and the results obtained 
were identical (up to four decimals) to the ones in the n = 100 column. 



In particular X is an additive process with representation ()4.4|) given by /? = 0, a = and the 
density of the compensator given by 

7 / -h'{x/t^), x>0, 

oft, X) = —r— < 

\ h'ix/t-y), x<o. 

In the specific case of the VG-Sato process we take the function h in ()7.14p to equal 

(7.15) h{x) := C (exp(-G|x|)I{^<o} + exp(-Mx)I{^>o}) , where M > 2T^ , C,G>0 

and T is the maturity of interest. Note that the analogue of the integral in (|5.11|) in the current 
setting can be computed in terms of the function h as: 

g{t,y)dy = j{h{u/r)-h{v/t^)) for 0<u<v. 

An analogous expression holds u < v < 0. Numerical results are contained in Table [721 

8. Conclusion 

In this paper we presented an algorithm for pricing barrier options in one-dimensional Markovian 
models based on an approximation by continuous-time Markov chains. The approximate barrier 
option prices are obtained by calculating the corresponding first-passage distributions for this chain. 
To illustrate the flexibility of the method we implemented the algorithm for a number of diffusion 
and exponential Levy models, a local volatility model with jumps and a model with time-dependent 
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VCr-Srito nrorpss 


n = 50 


n = 500 


n = 1000 


L\ — IDUU 


u.zoui 


n 9'?n/i 

U.ZoU4 


n 9'?n/i 


Ly — ZUUU 


U./oUZ 


U.ZoUO 


U.ZoUO 






U.ZoUD 


n 9'^nfi 

U.ZoUO 


N=2800 


0.2304 


0.2307 


0.2307 


N=3200 


0.2304 


0.2307 


0.2307 



Table 7. The parameter values for the VG-Sato process are taken from ^10 where the model 
was calibrated to options on the Amazon stock. The values are v = 0.7077, 7 = 0.4465, = 
— 1.13540, a = 0.7721 and the parameters C,G,AI in (|7.15p are given by the formulas 1/C = v, 
1/G = (-^6*^2/4 + cr2iy/2 - 9i^/2), 1/M = + + eu/2). The market data and 

contract details for the double-no-touch option are So = 100, r — 0.02, d — 0, T — 1/12 and 
i = 80, u = 120. The state-space of the approximating chain contains A'^ points and n denotes the 
number of time-steps. We also computed the double-no-touch prices in the case n — 2000 for all 
N reported in this table and found that the numbers are the same (the first four decimals) as the 
ones in the column n = 1000. 

jump-distributions. In the cases of the diffusion and Levy models, where results had been obtained 
before in the literature, the algorithm produced outcomes that accurately matched those results, 
while in the other cases numerical convergence was shown. We also provided a mathematical proof 
of the convergence of this algorithm to the true prices. However, to assess a priori the accuracy 
of the results produced by the algorithm it would be required to establish error estimates and 
rates of convergence for this Markov chain approximation method, which is left for future research. 
Although in principle the method also applies to higher-dimensional Markov processes, the size of 
the generator matrix would make straightforward application of the algorithm too time-consuming. 
Investigation of this extension is another topic left for future research. 

Appendix A. Proofs 

A.l. Proof of Lemma [Tl Without loss of generality we can restrict to the case r = 0. To prove 
the assertion we need to show that, for / G I? C Co(E), gt{x) ^ as t | 0, where 

gt{x) := t-\E,[f{St)I{t<r^y]+E,[f{Sr^)^t^,^y]-f{x))-k{x), 

and k is defined in (j2.12p . By definition, gt{x) = for x £ A whereas for x G ]E\^, 

gt{x) = t-HE,[f{St)l{t<r^}]+E,[f{SrJI{t>r^y]-f{x))-Cf{x) 
(A.l) = {t~\BM{St)] - fix)) - Cf{x)} - {t-iE,.[(/(5i) - fiS^J) Iit>r^}]} , 

which tends to zero as t | 0. Indeed, note that the first term in (IA.1|1 tends to zero since 

(A.2) Ct := sup \t~\-E,[f (St)] - fix))- £f{x)\^0 as U 0. 

x&E 
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Furthermore, the strong Markov property and the fact that Pq/ = / imply that the second term 
tends to zero, as follows: 

\t~'E, [I{i>,^} {Pt-rJ (5rJ - / (SrJ)] \ < P. (t^ < t) SUp S-'\P.J{x) - f{x)\. 

s<t,x£E 

The latter tends to zero as t | 0, since the second term is bounded in view of (IA.2p . 

Note that if / is such that liuix^dA ^f{x) = 0, then k G Co(IE). For such an /, if is a 
Feller process, the point-wise convergence for every x G E of t~^{P^f{x) — f{x)) to k{x) as t | 
implies that C^f = k, where is the generator of P^. This fact follows as a consequence of the 
Hille-Yosida theorem, see e.g. [Ml Lemma 31.7]. 



A. 2. Proof of the Feller property and the form of the generator of an exponential 
additive process. Let S be an exponential additive process as given in (I4.3p ~ ()4.4p and let (D, Y) 
be as in ()4.ip . In view of the right continuity of t {Dt,St), the boundedness of the function 
{s,x) ^ f{s,x) and the dominated convergence theorem, it follows that Qtf{s,x) converges to 
f{s, x) for every s £ [0, T] and x € E as t | 0. The linear function g Qtg maps Co([0, T] x E) to 
itself, which follows by combining the dominated convergence theorem, the spatial homogeneity of 
X and the fact that g £ Co([0,T] x E). Thus, {D,S) is a Feller process. 

Denote the marginal distribution and corresponding characteristic function of the increment 
Xt-i-s — Xg by fJ-s,s+t and fls^s+t{z) respectively. In view of the form (j4.4p and the independent 
increments property of X it holds that 




ipu{z)du where 



z f 

^puiz) := izf3{u) - —a{uf + / [e'^^ - 1 - iyz]g{u, y)dy. 

The characteristic function of a compound Poisson process with Levy measure t~^fj,s^s+tn is 
given by 

Jln{z) = expl^ {e''='-l)fis,s+tMx)\ 

\ n J M. / 

= exp{t-^{-j2s^s+tJz) -1)) 

= exp(t;:i(e/^'"^"(^)'i--l)). 

For fixed z and any sequence (tn)nGN that converges monotonically to zero, we find that, as n ^ oo, 

exp{ips{z)) =: ^s{z), 

since t^j^ il) s+u{z)d.u is right differentiable in zero with derivative '4>s{z)i in view of the dominated 
convergence, the integrability condition (14. 5p and the continuity of (3, a and g{-,y). Here Jlsiz) 
denotes a characteristic function of an infinitely divisible distribution. Hence an argument analogous 
to that in the proof of Theorem 31.5 in Sato [H] can be constructed to verify that for / G C|(M-|-) 
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the identity holds 

2 / \ 2 

limt-\Qtf-f)is,x) = p{s)xVfix) + ^^^Afix) 

+ [ [fixe') - fix) - I{|,|<i|V/(x)(e^ - l)]gis, z)dz, 
Jr 

where the limit is taken in the strong sense. 

A. 3. Proof of convergence. In this proof we shall employ standard convergence results for 
Markov processes that can be found in Ethier & Kurtz |19|. 

The proof is based on the following results. Let Ptfix) = Ea;[/(Xt)] and Pj;^^ fnix) = E^[/„(Xj^"-')] 
where, for any function / : E ^ M, we write /„ to denote /„ = /|g„- Denote by 2?(E) a core of the 
infinitesimal generator C, and write gn ^ g sup^g^j^ Idnix) — gix)\ ^ as n — > oo. Then it holds 
that (Ethier and Kurtz [lUt, Theorem 1.6.1]) 

A^"^/n^/:/ forall/eP(E) 

implies that 

(A.3) p/"Vn -^Ptf for ah / G Co(E) and t > 0. 

Furthermore, if (jA.Sp holds, then for any starting point Sq = Xq^^ = x E E, n G N, we have 
S, that is, converges weakly to S on -D(M), the Skorokhod space of cadlag real- valued 
functions endowed with the Skorokhod topology (Ethier and Kurtz [191 Theorem 4.2.11]). 

Thus, condition (16. 2p implies that X^") S, or, equivalently, for any continuous bounded 
function 5 : i:>(M) ^ R it holds that 

(A.4) E,[<7(X("))] ^ E,[giS)] as n ^ 00. 

In fact (|A.4p holds if g : D(]R) ^ M is bounded and continuous on some subset C of Z)(]R) that 
satisfies Px[S £ C] = 1 (see Jacod and Shiryaev [29l Section VI.3a]). 

To complete the proof we thus have to establish the continuity of the barrier payoff, which is for 
any uj £ Z)(M) given by 

a; ^ giuj) = e-'-(^^(-)^^)/(u;(r^(a;) A T)), 

in the Skorokhod topology. We refer to Jacod and Shiryaev [29] for background on the Skorokhod 
topology. Define 

Ti^uioj) := inf{t > : u;(t) or w(t-) ^ (^,u)}, 
T^^iu) := mi{t>0:uit)(^[i,u]}, 

Jiuj) := {s G M+ : a;(s) / a;(s-)}, 

Viiv) := {y = (yi,y2):i;(^)<r+(a;)}, 

V'iu) := {y = (yi,y2) : ^^(c^) G J(a;), u^iTyiu)-) G {2/1,2/2}}. 
Then the following general result holds true: 
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Lemma 2. (a) At every uj such that t ^ -/(w) the map D(M.) M? given by 

u! I— > I inf oj{s) A 0, sup oj{s) V ) 
\0<s<t 0<s<t J 

is continuous. 

(h) At every uj such that {i,u) ^ V{io) and T^^u(a;) < oo the map D{M) M-|_ given by 
is continuous. 

(c) At every to such that t ^ J{^) and {£,u) ^ V{uj) U V'{uj) the map D(R) M given by 

to 1-^ uj{T£^ui'-^) A t) 

is continuous. 

The proofs of Lemma 2(a), (b) and (c) are straightforward adaptations of Propositions VL2.4, 
VI.2.10 and VL2.11 in Jacod and Shiryaev [29]. 

Assume now that the process S" is a coordinate process on the canonical probabihty space -D(M), 
i.e. UJ = S{uj) for each uj G D(R). The fact that S is quasi-left continuous (as it is a Feller process, 
e.g. Thorem 2.4]) implies that at each fixed time t the path t St{uj) is continuous almost 
surely, that is, Px{t ^ Ji^)) = 1 for any x Furthermore, quasi-left continuity also implies that 
Px[(yi,y2) G V'{uj)] = Px[Ty{uj) £ J{uj), uj{Ty{uj)-) G {^1,^2}] = for any pair yi,y2- Indeed, 
on the event {uj{Ty{uj) — ) G {yi,y2},Ty{u>) < 00}, it holds that uj{Ty{uj)—) = uj{Ty{uj)) G {yi,y2} 
almost surely, as for any increasing sequence of stopping times r„ < Ty converging to Ty it holds 
that uj{Tn) — > uj(Ty—) = uj{Ty) almost surely on {Ty{uj) < 00}. Thus, almost surely T^^ is equal to 

T>,„(u;) := inf{t > : w(t) i {l,u)}. 

Further, the condition (j6.ip implies that 

Px(rA = r) = 0. 

Indeed, the event {ta = T} is equal to the union of the events {St = Sr^ G {i,u}}, and {St = 
St-^ G (0,^) U {u,oo)}, both of which have zero probability; the former by condition (16. ip and the 
latter since Px{St / 5't-) = by quasi-left continuity of 5. 

Moreover, note that Assumption [2] implies that Pxi{£,u) ^ V{uj)) = 0. The proof of Theorem [2] 
is then completed by combining the foregoing with the convergence in equation ()A.4p and Lemma[2l 

A. 4. Proof of time-dependent convergence. Theorem [2] remains valid for a higher-dimensional 
Feller process if the first passage is defined on one of the coordinates, as the above proof carries 
over. The question we need to answer is how to construct the sequence of Markov chains that 
approximates the two-dimensional Feller process {D, Y) (see (14. ID for definition) in such a way that 
the limits in Corollary [2] hold. 
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A. 4.1. Construction of an approximating two-dimensional time-homogeneous chain. To fit the case 
of a time-inhomogeneous Markov process into this framework we will consider an approximation of 
the time-space process {D,Y). We approximate the time-space Feller Markov process {D,Y) by a 
two-dimensional (time-homogeneous) Markov chain (Z, X) living on a time-space grid T x G where 
T = {0} U with 6m = T/ {nm) and n, m G N. To define this chain, let C be the infinitesimal 
generator of {D,Y) and Ctf the corresponding space restricted generator. Assume that A^*\ 
i = 1, ... ,71 are generator matrices approximating the generator CiTin- Then the generator of the 
two-dimensional chain (Z, X) is specified as 



(A.5) A("'-)/(f, X) = 6-\fit + 5m, x) - fit, x)) + ^(A»/t)(x)I|(,„i 



)m5,„<t<im5m} 



i=l 



for (t, x) € T X G, where ft{x) = f{t, x), for any function / : T x G ^ M. 

The corresponding stochastic dynamics are described as follows. The first component is given 
by Zt = 6m ■ Nt, 6m times a standard Poisson process N with rate A = l/6m- Thus, Z is piecewise 
constant and moves by positive jumps of size 6m that occur after independent exponential times with 
mean A. For any t > 0, Zt follows a Poisson distribution with mean t and variance 6mt = tT/{nm). 
Hence, as n — > oo, Zt tends to t for every fixed t > 0. In fact, since {Zt — t}t>o is a martingale, 
Doob's maximal inequality implies that 



E 



sup {Zt - ty 

t<T 



< YaTlZx] = T'^/{nm), 



so that {Zt}t^[Q^T] converges uniformly to the deterministic unit drift {Dt}te[o,T]j almost surely. 
Conditional on Zt taking a value in [(i — l)6m,i6m), X evolves as a Markov chain with state-space 
G and generator matrix A*^*). 

In this setting we consider the barrier option with a randomised time of maturity 

rpinm) ^ -j^fi^ >0: Zt = T}, 

the first time that Z hits the level T = mn6m, which is equal to the nm-th jump time of Z. 
Thus T("™) is distributed as the sum of mn independent exponential random variables with mean 
T /{nm), and follows a Gamma distribution with mean T and variance T^/ {nm). The corresponding 
(stochastic) discounting is defined as 

(A.6) Rt= I p{Zs)<^s, where /5(^i) = V I{(i-i)m5„<«<im5,„}- 

•^0 i=l 

The idea of randomising the maturity was employed before by Carr [7] to approximate a finite 
maturity American put option by replacing its maturity by an independent Gamma(n, T/n) random 
variable, and this technique has also been employed to value barrier- type options (see e.g. [5]). 
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Corollary 3. For any T > 0, p in (1A.6P and any functions (/> : G — > M and -0 : G — > M with 
ip{x) = for X £ G, it holds that 



(A.7)E,. 0(X 



(A.8)E, 



'J'(nm) ) -'-|T->J'(nm)| 



^(A,)I 



/ - 5™, a(i) 



J/-5^A(") 



X G 



a; G G, 



where 6m = T/ {nm) and the stopping time t is defined in ()3.3p . 

Proof. The strong Markov property of (Z, X) at the stopping time T^^^ yields that 



E, 



E, 



I{r>r(i)}Ex^(i) 



{r>r('""-i)} 



2^(mn — 1) 



smce is in distribution equal to the sum of i independent exponential random variables. It 
follows from Theorem [1] that for any / : G ^ R 

POO 

Ex[W(i)}/(^T(i))] = / Xe-^'ES{r>t}f{Xt)]dt 

Jo 

= [(/- A-^AW)^V](a^), where A 



The expression (IA.7P then follows by induction. The proof of the identity (lA.Sp is similar and is 
therefore omitted. |--| 

A. 4. 2. Proof of Corollary\^ For any sequence {m{n))n&% of natural numbers, let (Z^™''"^ 
be a sequence of continuous-time Markov chains with state-spaces T^"") xG^"'^ and generator matrices 
A^"'™^ as constructed in the previous section, such that their generators uniformly converge to the 
two-dimensional generator C of {D,Y), in the sense that 

A("'"^)/„(t,x)-£7(i,x) 



(A.9) 



(^n,m{f) = max 







for all functions / in a core of where /„ = /|[ot]xG(")- view of the form of the generator 
(lA.Sp and the condition (16. 3p in Corollary [21 it follows that (lA.Op holds true. If S satisfies (16. ID 
and Assumption [2] (with replaced by Po,x)! then the conclusions of Theorem [2] also apply to the 
two-dimensional Feller process {D,Y) and the sequence (Z^"*'"), (since the proof remains 

valid for this case). As a consequence the expectations (jA.7p - (|A.8p converge to the barrier option 
prices under the limiting model (D,Y). To complete the proof we must check that the sequences 
(lAr7ll - (fX8]l converge to the same hmits as (fO]) ~ (fiT0]l . 

The latter follows by observing that, for a given m G N, a matrix A, a maturity T > and 
6m = T/{nm) it holds that, 



exp < —A 
[ n 

{i-6m^r 



n 2n^ Zw^m 



A' 



on 



as n tends to infinity. Hence, for given matrices , we have that 



(A.IO) 



exp 



.Ad) 



n 



■ ■ ■ exp 



n 



A(") - (I - 5™A(i))— •••(/- (^^AW)- 



< 
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where C*-"-* is some constant that depends on A'^^^ . . . , A*-") and || • || is a matrix norm. 

Take the sequence of natural numbers (m(n))„gp^ such that it satisfies m{n) > C(") for all n e N, 
where {C^^^)neN sue the constants in (|A.10p . that correspond to the sequence {Z^"^'"'\ X^"'^). It 
then follows from (TOOl) that the expectations (IXTjl - dXSll and (liTOll- diTOl) converge to the same 
limit as n ^ oo. 

Appendix B. Non-uniform state-space of the Markov chain X 

The algorithm for constructing a non-uniform grid {xi, . . . , xn} with G 2N points for the 
triplet a < s < b with density parameters gi and g2 is given by the following procedure. 

(1) Compute ci = arcsinh ^^^)) C2 = arcsinh (^^^ 

(2) Define the lower part of the grid by the formula Xk ■= s + 5isinh(ci(l — (A; — 1) /{N/2 — 1))) 
for k G {1, . . . , A^/2}. Note that xi = a, xj\f/2 = •S- 

(3) Define the upper part of the grid using the formula Xi^^pf/2 '■= s + g2smh{c22k / N) for 
k € {1,..., N/2}. Note that xn = u. 

The non-uniform state-space G C E for the Markov chain X, used for the pricing of a double- 
barrier option with barrier levels I < u, can now be constructed in the following way. Pick A^ G 3N, 
choose the smallest xi and largest xn values of the grid G and apply the procedure above three 
times to the triplets xi < I < {s — l)/2, {s — l)/2<s< {u — s)/2 and {u — s) /2 < u < xjy, each time 
with N/3 points and some density parameters. The state-space G is obtained by concatenating the 
three grids (see also the Matlab function that generates this grid in [43j ) . 
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